rtsVis is a package with a single purpose: The creation of animated visualizations of time satellite imagery time series.
In this overview, we show how to create such visualizations on the example of the MODIS MCD43A4 product of the Western Cape, South Africa.
To get started, we load a number of packages. To begin with, we need rtsVis itself. We use raster to load and plot rasters and sf for vector objects.
rtsVis also goes hand in hand with moveVis from which we use some core functions. We also use magrittr for the convenient pipes.
library(rtsVis)
library(raster)
library(sf)
library(moveVis)
library(magrittr)
Perhaps most importantly, we need the data itself: A time series of raster images.
rtsVis works on lists of objects - lists of raster objects or lists of frames (plots). Most functions of the package take a list as input, modify it, and return the modified list.
The first step in a rtsVis pipeline is to get our time series into this format.
We can load the images using the raster::stack function.
modis_filenames <- list.files("Data/Modis_images/",pattern = ".tif$",full.names = T)
r_list <- lapply(modis_filenames,stack)
Or maybe we have a time series prepared in some other format. Sometimes, you will already have done some preprocessing, such as resampling or masking. What matters is that the images are matching in projection and extent. We have such a time series already stored as .rds file.
r_list <- readRDS("Beispieldaten/MODIS/WesternCape/MODIS_WesternCape.rds")
First, we take a look at one of our images.
plotRGB(r_list[[20]],1,4,3,stretch="lin")
We can see that it looks quite good. But there are very small patches of what appear to be missing values.
To fix them, we use the first rtsVis function: ts_fill_na(). This uses the raster::approxNA() function to fill NA values from the closest temporal neighbor in the time series - a crude approach, but reasonably quick.
r_list_filled <- rtsVis::ts_fill_na(r_list)
plotRGB(r_list_filled[[20]],1,4,3,stretch="lin")
Much better! In the southwest, out on the Atlantic Ocean, we have an area which cannot be filled because we have no values for it in any of our MODIS images.
Our MODIS Time Series is already at a very good temporal resolution of one day. Great for analysis - but for visualizations, we do not want to be restricted by the temporal resolution of our source.
We will use the rtsVis::ts_raster function to linearly interpolate our time series to a higher temporal resolution. For this, we need:
r_listr_times for each input imageout_times for each desired output imageOur input images we have already. So all we need are the times.
How we get the times depends on the source of our images. Here, we stored our acquisition dates in the names of our list object. In other cases, we might read them out of the metadata or the file names. However we get there, what matters is that we get a vector of POSIX objects of the same length as our list of images, containing the input time for each image.
head(names(r_list_filled))
## [1] "2019-07-01" "2019-07-02" "2019-07-03" "2019-07-04" "2019-07-05"
## [6] "2019-07-06"
in_dates <- as.POSIXct(names(r_list_filled))
head(in_dates)
## [1] "2019-07-01 CEST" "2019-07-02 CEST" "2019-07-03 CEST" "2019-07-04 CEST"
## [5] "2019-07-05 CEST" "2019-07-06 CEST"
Our output times are all the dates for which we will receive an interpolated output image. The output times, therefore, should also be a vector of POSIX (although they need not match the input images). How we get them is, again, up to us. If we only want to increase the temporal resolution, we can simply create a sequence of POSIX, like this:
out_dates <-seq.POSIXt(from = in_dates[1],
to = in_dates[length(in_dates)],
length.out = length(in_dates)*2 )
head(out_dates)
## [1] "2019-07-01 00:00:00 CEST" "2019-07-01 11:57:51 CEST"
## [3] "2019-07-01 23:55:42 CEST" "2019-07-02 11:53:33 CEST"
## [5] "2019-07-02 23:51:24 CEST" "2019-07-03 11:49:15 CEST"
We pass our images and times to ts_raster. The option fade_raster activates linear interpolation instead of nearest-neighbor interpolation.
r_list_filled_interpolated <- rtsVis::ts_raster(r_list = r_list_filled,
r_times = in_dates,
out_times = out_dates,
fade_raster = T)
length(r_list_filled)
## [1] 154
length(r_list_filled_interpolated)
## [1] 308
We have doubled the length of our time series. We should bear in mind that the interpolated values no longer correspond to actual measurements and should not be used for analysis. Also, any further processing will take longer. So we should use this power responsibly. For visualizations, however, this is a great boon, as we have doubled the frequency of our images.
Also, the resulting list has the time written in its attributes. This is crucial for the application of some rtsVis functions.
head(rtsVis:::.ts_get_frametimes(r_list_filled_interpolated))
## [1] "2019-07-01 00:00:00 CEST" "2019-07-01 11:57:51 CEST"
## [3] "2019-07-01 23:55:42 CEST" "2019-07-02 11:53:33 CEST"
## [5] "2019-07-02 23:51:24 CEST" "2019-07-03 11:49:15 CEST"
Now we have a time series in our required format and at our desired resolution. We can go on to create some visualization.
All the visualizations we create will be ggplots. From our list of raster objects will come a list of ggplot frames, one for each input raster object.
ts_makeframes(). These can be further enhanced by adding
moveVis::add_northarrow()rtsVis::ts_add_positions_to_frames()The returned list of frames can be animated using moveVis::animate_frames().
Flow frames visualize amounts, distributions or proportions of raster values within an area. They can, in principle, take many forms such as lines, histograms, boxplots, scatterplots, pie charts or barcharts. rtsVis currently implements violin plots and line plots. However, custom plot types can be used as well.
By default, these plots will be calculated for the entire area covered by the raster. However, if we want to distinguish between different areas of interest within our study area. In this case, we can provide the spatial objects and the plots will use the spatial objects as a grouping variable.
Here, we use four administrative areas which lie within out study area.
polygons <- st_read("Beispieldaten/Ancillary/WesternCape/Western_Cape_polygons.shp") %>% st_transform(crs = st_crs(r_list_filled_interpolated[[1]]))
#Plot the polygons
plotRGB(r_list_filled_interpolated[[50]],1,4,3,stretch="lin")
plot(st_geometry(polygons),add=T,col="red",alpha=0.5)
We can see that they are of different sizes and shapes and cover different land uses. This may be important when evaluating our plots later. But first, we create them. Only the first three arguments are truly necessary, however, we make use of the option to customize the plots.
violin_frames <- rtsVis::ts_flow_frames(r_list = r_list_filled_interpolated,
positions = polygons,
plot_function = "violin",
# Optional Arguments
position_names = polygons$NAME_3,
band_names = c("620 - 670","841 - 876","459 - 479","545 - 565"),
band_colors = c("firebrick3","darkorchid3","dodgerblue3","olivedrab3"),
band_legend_title = "Wavelength [nm]",
position_legend = F,
legend_position = "left",
band_legend=F,
aes_by_pos = F)
The result is a list of the same length as our list of input rasters.
length(violin_frames)
## [1] 308
And each object in our new list is a ggplot. We can take a look at an individual violin_frame.
violin_frames[[50]]
## Warning: Removed 68 rows containing non-finite values (stat_ydensity).
In the same way we used ts_flow_frames() to create a violin plot of polygon data, we can use it to create a line plot of data under a set of points.
points <- st_read("Beispieldaten/Ancillary/WesternCape/Western_Cape_pts.shp") %>% st_transform(crs = st_crs(r_list_filled_interpolated[[1]]))
line_frames <- rtsVis::ts_flow_frames(r_list = r_list_filled_interpolated,
positions = points[c(1,3,4),],
plot_function = "line",
# Optional Arguments
position_names = points$pointname[c(1,3,4)],
band_names = c("620 - 670","841 - 876","459 - 479","545 - 565"),
band_colors = c("firebrick3","darkorchid3","dodgerblue3","olivedrab3"),
band_legend_title = "Wavelength [nm]",
position_legend = T,
legend_position = "left",
band_legend=T,
aes_by_pos = T)
line_frames[[length(line_frames)]]
## Reading layer `Western_Cape_pts' from data source `D:\DEV\TSVis\Beispieldaten\Ancillary\WesternCape\Western_Cape_pts.shp' using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 2 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 1682063 ymin: -3736239 xmax: 1797136 ymax: -3663320
## projected CRS: Sinusoidal
We could, in principle, also create a violin plot of the data extracted under points. Then we would probably want to use a buffer to our point data, like so:
points <- st_read("Beispieldaten/Ancillary/WesternCape/Western_Cape_pts.shp") %>% st_transform(crs = st_crs(r_list_filled_interpolated[[1]]))
line_frames <- rtsVis::ts_flow_frames(r_list = r_list_filled_interpolated,
positions = points[c(1,3,4),],
plot_function = "violin",
pbuffer = 1337, ##Use values in a circular buffer of 1337
##map units around the points
#
position_names = points$pointname[c(1,3,4)],
band_names = c("620 - 670","841 - 876","459 - 479","545 - 565"),
band_colors = c("firebrick3","darkorchid3","dodgerblue3","olivedrab3"),
band_legend_title = "Wavelength [nm]",
position_legend = T,
legend_position = "left",
band_legend=F,
aes_by_pos = T)
line_frames[[length(line_frames)]]
Violin and line charts are nice for some applications but often other plot types may be more suitable. Besides the options “line” and “violin”, we can also provide ts_flow_frames() with a custom plotting function.
There are a few requirements for such a function. It should:
Below, we create a plotting function that creates a barchart.
library(ggplot2)
barplot <- function(i, edf , pl,lp, bl, blt,plt, ps, vs,abp){
#The data up to the current frame (this will be plotted)
x = edf[edf$frame == i,]
#All data (this sets the frame)
y=edf
## generate base plot
p <- ggplot(x, aes(x = band,y=value,fill=band))
## style it
p <- p +
geom_bar(stat = "identity")+
theme_bw() +
theme(aspect.ratio = 1) +
scale_y_continuous( breaks = vs,limits = c(min(vs),max(vs)))+
theme(legend.position = lp)+
facet_grid(~position_name)+
scale_x_discrete(guide = guide_axis(n.dodge = 2))
#add the colors
p <- p +
scale_fill_manual(values = x$band_colors,breaks = x$band, name=blt)
return(p)
}
We can then pass this function to ts_flow_frames().
bar_frames <- rtsVis::ts_flow_frames(r_list = r_list_filled_interpolated,
positions = points[c(1,3,4),],
plot_function = barplot,
position_names = points$pointname[c(1,3,4)],
band_names = c("620 - 670","841 - 876","459 - 479","545 - 565"),
band_colors = c("firebrick3","darkorchid3","dodgerblue3","olivedrab3"),
band_legend_title = "Wavelength [nm]",
position_legend = T,
legend_position = "left",
band_legend=T,
aes_by_pos = T)
bar_frames[[50]]
Many other plot types are possible. For inspiration, visit the r-graph-gallery!
Also visit Claus Wilke’s excellent book which deals with the most common pitfalls of data visualization.
Spatial frames are, quite simply, images that can be discrete, gradient or RGB, depending of the type of input raster.
The frames can be created using the ts_makeframes function.
Two things to note here:
r_frames <- rtsVis::ts_makeframes(x_list = r_list_filled_interpolated,
l_indices = c(1,4,3),
minq = 0.2,
maxq = 0.999)
## [1] "Guessing raster type."
## [1] "Detected 3+ layers, choosing raster type 'RGB'"
r_frames[[50]] # Examine a single frame
We can see that our frames are plotted as ggplots with latitude and longitude as the axes.
We can use functionality from the moveVis package to enhance our images with additional elements.
r_frames_styled <- r_frames %>%
moveVis::add_labels(x = "Longitude", y = "Latitude")%>%
moveVis::add_northarrow(colour = "white", position = "bottomright") %>%
moveVis::add_timestamps(type = "label") %>%
moveVis::add_progress()
r_frames_styled[[50]] # Examine a single frame
Now the spatial and temporal context of the map is greatly increased. In some cases, we want to go even further.
We may want to visualize positions, certain points or areas of significance on the map. This can be useful to visualize the areas of interest we explore in the chats. In our case here, we use the same positions we previously used to create the flow frames.
In our first case, we add the municipalities - which we used to create the violin plots - and add them as polygons to our frames.
r_frames_with_polygons <-
rtsVis::ts_add_positions_to_frames(r_frame_list = r_frames_styled,
positions = polygons,
aes_by_pos = T,
legend_position = "right",
position_names = polygons$NAME_3,
position_legend_title = "Municipality",
psize = 1,
add_text = T,
tsize = 3,
t_hjust = 5000,
ttype = "label")
r_frames_with_polygons[[50]] # Examine a single frame
In our second case, we do it with our previously used points.
r_frames_with_points <-
rtsVis::ts_add_positions_to_frames(r_frame_list = r_frames_styled,
position_names = points$pointname,
positions = points,
add_text = T,
tsize = 4,
ttype = "label",
t_vjust = 8500)
r_frames_with_points[[50]] # Examine a single frame
we can use moveVis::join_frames to join frames. This way we can combine our charts which show the distribution of values within our areas of interest with a spatial element that presents their spatial context. See the join_frames() documentation for more details.
point_joined<- moveVis::join_frames(frames_lists = list(r_frames_with_points,line_frames))
point_joined[[10]]
Lists of frames, as returned by ts_makeframes(), ts_flow_frames(), or moveVis::join_frames(), can be animated in a number of different formats ùsingmoveVis function, moveVis::animate_frames(). See the documentation of the function for more details. The resulting file is written to disk.
moveVis::animate_frames(point_joined,
out_file = "WesternCape_MODIS_point.gif",
width=1600)
To learn more about creating custom plot types, access the guide and test material at the associated github repo.
For creating visualization with movement data, visit moveVis.
For inspiration, visit the r-graph-gallery!
To learn about data visualization see Claus Wilke’s excellent book.